
## @knitr AALEN_SETUP

rm(list=ls())

#library(foreign)
#library(stringr)
#library(xlsx)

# from \url{https://www.r-bloggers.com/identifying-the-os-from-r/}
get_os <- function(){
  sysinf <- Sys.info()
  if (!is.null(sysinf)){
    os <- sysinf['sysname']
    if (os == 'Darwin')
      os <- "osx"
  } else { ## mystery machine
    os <- .Platform$OS.type
    if (grepl("^darwin", R.version$os))
      os <- "osx"
    if (grepl("linux-gnu", R.version$os))
      os <- "linux"
  }
  tolower(os)
}

get_machine  <- function(){
    sysinf  <- Sys.info()
    if (!is.null(sysinf)){
        my.machine  <- sysinf['nodename']
    }
    tolower(my.machine)
}


if(get_os()=="windows" & get_machine()=="mompracen") {
    path <- setwd('C:/Users/giaco/Documents/MEGAsync/')
} else if (get_os()=="windows" & get_machine()!="mompracen") {
    path <- setwd('C:/Users/grieco/Dropbox/')
} else if (get_machine()=="aramis") {
    path <- setwd('/home/giacomo68/Dropbox/')
} else {
    path <- setwd('/home/giacomo/Dropbox/')
}


library(survival)
library(xtable)
library(ggplot2)
library(english)
library(Publish)
library(gridExtra)
library(ggpubr)
library(ggfortify)
library(mice)
library(timereg)
library(reshape2)
library(pammtools)
library(patchwork)

Data <- read.csv('gc-jg-project/kampala/data/data-kampala.csv',
                 stringsAsFactors=FALSE)



Data$country_id <- as.numeric(factor(Data$country_name))

m <- 10

imp <- mice(Data[,c("kampala",
                 "sidea.1991",
                 "sideb.1991",
                 "my.mil.cap.lag",
                 "v2x_polyarchy.lag",
                 "v2x_libdem.lag",     
                 "v2x_partipdem.lag",
                 "v2x_delibdem.lag", 
                 "v2x_egaldem.lag",
                 "gdp.cap.lag",
                 "pop.lag",
                 "common.law",
                 "year",
                 "rule.law.lag",
                 "aid.recipient",
                 "aid.recipient.q",
                 "All_NGO.lag",
                 "country_id")], m = m, maxit = 10, seed = 1234)

models <- list()
model_results <- list()
merged_imputations <- list()
m1.model_results <- list()
m2.model_results <- list()
m3.model_results <- list()
m4.model_results <- list()
m5.model_results <- list()
m6.model_results <- list()


# Initialize lists to store all models
all_models <- list(
  m1 = vector("list", m),
  m2 = vector("list", m),
  m3 = vector("list", m),
  m4 = vector("list", m),
  m5 = vector("list", m),
  m6 = vector("list", m)
)


for (i in 1:m) {
    complete_data <- complete(imp, i)
    complete_data$.imp <- i
    
    # Merge with original data
    merged_data <- merge(complete_data, 
                       Data[,c("country_id","year","country_name",
                              "kampala.t0","kampala.t1",
                              "sidea.1945","sideb.1945",
                              "sidea.force.1945","sideb.force.1945",
                              "sidea.no.force.1945","sideb.no.force.1945",
                              "sidea.force.1991","sideb.force.1991",
                              "sidea.no.force.1991","sideb.no.force.1991")], 
                       by = c("country_id", "year"), 
                       all.x = TRUE)
    
    # Create transformed variables
    merged_data$"Common law" <- merged_data$common.law
    merged_data$"Polyarchy" <- merged_data$v2x_libdem.lag
    merged_data$"GDP capita (log)" <- log(merged_data$gdp.cap.lag)
    merged_data$"Population (log)" <- log(merged_data$pop.lag)
    merged_data$"Milex per cap" <- merged_data$my.mil.cap.lag

    
    # Model 1: 1991 MIDs (all)
    merged_data$"MID side A" <- merged_data$sidea.1991
    merged_data$"MID side B" <- merged_data$sideb.1991
    all_models$m1[[i]] <- aalen(Surv(kampala.t0, kampala.t1, kampala) ~
                               `MID side A` + `MID side B` + `Milex per cap` +
                               `Polyarchy` + `GDP capita (log)` + 
                               `Population (log)` + `Common law`,
                               data = merged_data, 
                               clusters = Data$country_name,
                               id = Data$country_name,
                               start.time = min(Data$kampala.t0, na.rm = TRUE))
    
    # Model 2: 1991 MIDs (force only)
    merged_data$"MID side A" <- merged_data$sidea.force.1991
    merged_data$"MID side B" <- merged_data$sideb.force.1991
    all_models$m2[[i]] <- aalen(Surv(kampala.t0, kampala.t1, kampala) ~
                               `MID side A` + `MID side B` + `Milex per cap` +
                               `Polyarchy` + `GDP capita (log)` + 
                               `Population (log)` + `Common law`,
                               data = merged_data, 
                               clusters = Data$country_name,
                               id = Data$country_name,
                               start.time = min(Data$kampala.t0, na.rm = TRUE))


    ## Model 3: 1991 MIDs (no force)
    merged_data$"MID side A" <- merged_data$sidea.no.force.1991
    merged_data$"MID side B" <- merged_data$sideb.no.force.1991
    all_models$m3[[i]] <- aalen(Surv(kampala.t0, kampala.t1, kampala) ~
                               `MID side A` + `MID side B` + `Milex per cap` +
                               `Polyarchy` + `GDP capita (log)` + 
                               `Population (log)` + `Common law`,
                               data = merged_data, 
                               clusters = Data$country_name,
                               id = Data$country_name,
                               start.time = min(Data$kampala.t0, na.rm = TRUE))


    # Model 4: 1945 MIDs (all)
    merged_data$"MID side A" <- merged_data$sidea.1945
    merged_data$"MID side B" <- merged_data$sideb.1945
    all_models$m4[[i]] <- aalen(Surv(kampala.t0, kampala.t1, kampala) ~
                               `MID side A` + `MID side B` + `Milex per cap` +
                               `Polyarchy` + `GDP capita (log)` + 
                               `Population (log)` + `Common law`,
                               data = merged_data, 
                               clusters = Data$country_name,
                               id = Data$country_name,
                               start.time = min(Data$kampala.t0, na.rm = TRUE))


    # Model 5: 1945 MIDs (force only)
    merged_data$"MID side A" <- merged_data$sidea.force.1945
    merged_data$"MID side B" <- merged_data$sideb.force.1945
    all_models$m5[[i]] <- aalen(Surv(kampala.t0, kampala.t1, kampala) ~
                               `MID side A` + `MID side B` + `Milex per cap` +
                               `Polyarchy` + `GDP capita (log)` + 
                               `Population (log)` + `Common law`,
                               data = merged_data, 
                               clusters = Data$country_name,
                               id = Data$country_name,
                               start.time = min(Data$kampala.t0, na.rm = TRUE))

    ## Model 6: 1945 MIDs (no force)
    merged_data$"MID side A" <- merged_data$sidea.no.force.1945
    merged_data$"MID side B" <- merged_data$sideb.no.force.1945
    all_models$m6[[i]] <- aalen(Surv(kampala.t0, kampala.t1, kampala) ~
                               `MID side A` + `MID side B` + `Milex per cap` +
                               `Polyarchy` + `GDP capita (log)` + 
                               `Population (log)` + `Common law`,
                               data = merged_data, 
                               clusters = Data$country_name,
                               id = Data$country_name,
                               start.time = min(Data$kampala.t0, na.rm = TRUE))

}


m1_coef <- lapply(all_models$m1, function(x) x$cum)
m1_coef_array <- simplify2array(m1_coef)
m1_coef_mean <- apply(m1_coef_array, c(1, 2), mean)

m1_vars <- lapply(all_models$m1, function(x) x$robvar.cum)
m1_vars_array <- simplify2array(m1_vars)
m1_vars_within <- apply(m1_vars_array, c(1, 2), mean)
m1_vars_between <- apply(m1_coef_array, c(1, 2), var)
m1_vars_total <- m1_vars_within + (1+1/m)*m1_vars_between

m1_coef_mean_long <- melt(m1_coef_mean, variable.name = "variable", value.name = "mean.coef")
m1_vars_total_long <- melt(m1_vars_total, variable.name = "variable", value.name = "total.var")
m1_aalen0 <- merge(m1_coef_mean_long,m1_vars_total_long,by=c("Var1","Var2"),all.x=TRUE)
m1_aalen0 <- m1_aalen0[order(m1_aalen0$Var2,m1_aalen0$Var1),]
m1_aalen_time <- m1_aalen0[m1_aalen0$Var2=="time",c("Var1","mean.coef")]
colnames(m1_aalen_time) <- c("Var1","Time")
m1_aalen1 <- m1_aalen0[m1_aalen0$Var2!="time",]
m1_aalen <- merge(m1_aalen1,m1_aalen_time,by="Var1",all.x=TRUE)
m1_aalen <- m1_aalen[order(m1_aalen$Var2,m1_aalen$Var1),]


m2_coef <- lapply(all_models$m2, function(x) x$cum)
m2_coef_array <- simplify2array(m2_coef)
m2_coef_mean <- apply(m2_coef_array, c(1, 2), mean)

m2_vars <- lapply(all_models$m2, function(x) x$robvar.cum)
m2_vars_array <- simplify2array(m2_vars)
m2_vars_within <- apply(m2_vars_array, c(1, 2), mean)
m2_vars_between <- apply(m2_coef_array, c(1, 2), var)
m2_vars_total <- m2_vars_within + (1+1/m)*m2_vars_between

m2_coef_mean_long <- melt(m2_coef_mean, variable.name = "variable", value.name = "mean.coef")
m2_vars_total_long <- melt(m2_vars_total, variable.name = "variable", value.name = "total.var")
m2_aalen0 <- merge(m2_coef_mean_long,m2_vars_total_long,by=c("Var1","Var2"),all.x=TRUE)
m2_aalen0 <- m2_aalen0[order(m2_aalen0$Var2,m2_aalen0$Var1),]
m2_aalen_time <- m2_aalen0[m2_aalen0$Var2=="time",c("Var1","mean.coef")]
colnames(m2_aalen_time) <- c("Var1","Time")
m2_aalen1 <- m2_aalen0[m2_aalen0$Var2!="time",]
m2_aalen <- merge(m2_aalen1,m2_aalen_time,by="Var1",all.x=TRUE)
m2_aalen <- m2_aalen[order(m2_aalen$Var2,m2_aalen$Var1),]


m3_coef <- lapply(all_models$m3, function(x) x$cum)
m3_coef_array <- simplify2array(m3_coef)
m3_coef_mean <- apply(m3_coef_array, c(1, 2), mean)

m3_vars <- lapply(all_models$m3, function(x) x$robvar.cum)
m3_vars_array <- simplify2array(m3_vars)
m3_vars_within <- apply(m3_vars_array, c(1, 2), mean)
m3_vars_between <- apply(m3_coef_array, c(1, 2), var)
m3_vars_total <- m3_vars_within + (1+1/m)*m3_vars_between

m3_coef_mean_long <- melt(m3_coef_mean, variable.name = "variable", value.name = "mean.coef")
m3_vars_total_long <- melt(m3_vars_total, variable.name = "variable", value.name = "total.var")
m3_aalen0 <- merge(m3_coef_mean_long,m3_vars_total_long,by=c("Var1","Var2"),all.x=TRUE)
m3_aalen0 <- m3_aalen0[order(m3_aalen0$Var2,m3_aalen0$Var1),]
m3_aalen_time <- m3_aalen0[m3_aalen0$Var2=="time",c("Var1","mean.coef")]
colnames(m3_aalen_time) <- c("Var1","Time")
m3_aalen1 <- m3_aalen0[m3_aalen0$Var2!="time",]
m3_aalen <- merge(m3_aalen1,m3_aalen_time,by="Var1",all.x=TRUE)
m3_aalen <- m3_aalen[order(m3_aalen$Var2,m3_aalen$Var1),]


m4_coef <- lapply(all_models$m4, function(x) x$cum)
m4_coef_array <- simplify2array(m4_coef)
m4_coef_mean <- apply(m4_coef_array, c(1, 2), mean)

m4_vars <- lapply(all_models$m4, function(x) x$robvar.cum)
m4_vars_array <- simplify2array(m4_vars)
m4_vars_within <- apply(m4_vars_array, c(1, 2), mean)
m4_vars_between <- apply(m4_coef_array, c(1, 2), var)
m4_vars_total <- m4_vars_within + (1+1/m)*m4_vars_between

m4_coef_mean_long <- melt(m4_coef_mean, variable.name = "variable", value.name = "mean.coef")
m4_vars_total_long <- melt(m4_vars_total, variable.name = "variable", value.name = "total.var")
m4_aalen0 <- merge(m4_coef_mean_long,m4_vars_total_long,by=c("Var1","Var2"),all.x=TRUE)
m4_aalen0 <- m4_aalen0[order(m4_aalen0$Var2,m4_aalen0$Var1),]
m4_aalen_time <- m4_aalen0[m4_aalen0$Var2=="time",c("Var1","mean.coef")]
colnames(m4_aalen_time) <- c("Var1","Time")
m4_aalen1 <- m4_aalen0[m4_aalen0$Var2!="time",]
m4_aalen <- merge(m4_aalen1,m4_aalen_time,by="Var1",all.x=TRUE)
m4_aalen <- m4_aalen[order(m4_aalen$Var2,m4_aalen$Var1),]


m5_coef <- lapply(all_models$m5, function(x) x$cum)
m5_coef_array <- simplify2array(m5_coef)
m5_coef_mean <- apply(m5_coef_array, c(1, 2), mean)

m5_vars <- lapply(all_models$m5, function(x) x$robvar.cum)
m5_vars_array <- simplify2array(m5_vars)
m5_vars_within <- apply(m5_vars_array, c(1, 2), mean)
m5_vars_between <- apply(m5_coef_array, c(1, 2), var)
m5_vars_total <- m5_vars_within + (1+1/m)*m5_vars_between

m5_coef_mean_long <- melt(m5_coef_mean, variable.name = "variable", value.name = "mean.coef")
m5_vars_total_long <- melt(m5_vars_total, variable.name = "variable", value.name = "total.var")
m5_aalen0 <- merge(m5_coef_mean_long,m5_vars_total_long,by=c("Var1","Var2"),all.x=TRUE)
m5_aalen0 <- m5_aalen0[order(m5_aalen0$Var2,m5_aalen0$Var1),]
m5_aalen_time <- m5_aalen0[m5_aalen0$Var2=="time",c("Var1","mean.coef")]
colnames(m5_aalen_time) <- c("Var1","Time")
m5_aalen1 <- m5_aalen0[m5_aalen0$Var2!="time",]
m5_aalen <- merge(m5_aalen1,m5_aalen_time,by="Var1",all.x=TRUE)
m5_aalen <- m5_aalen[order(m5_aalen$Var2,m5_aalen$Var1),]


m6_coef <- lapply(all_models$m6, function(x) x$cum)
m6_coef_array <- simplify2array(m6_coef)
m6_coef_mean <- apply(m6_coef_array, c(1, 2), mean)

m6_vars <- lapply(all_models$m6, function(x) x$robvar.cum)
m6_vars_array <- simplify2array(m6_vars)
m6_vars_within <- apply(m6_vars_array, c(1, 2), mean)
m6_vars_between <- apply(m6_coef_array, c(1, 2), var)
m6_vars_total <- m6_vars_within + (1+1/m)*m6_vars_between

m6_coef_mean_long <- melt(m6_coef_mean, variable.name = "variable", value.name = "mean.coef")
m6_vars_total_long <- melt(m6_vars_total, variable.name = "variable", value.name = "total.var")
m6_aalen0 <- merge(m6_coef_mean_long,m6_vars_total_long,by=c("Var1","Var2"),all.x=TRUE)
m6_aalen0 <- m6_aalen0[order(m6_aalen0$Var2,m6_aalen0$Var1),]
m6_aalen_time <- m6_aalen0[m6_aalen0$Var2=="time",c("Var1","mean.coef")]
colnames(m6_aalen_time) <- c("Var1","Time")
m6_aalen1 <- m6_aalen0[m6_aalen0$Var2!="time",]
m6_aalen <- merge(m6_aalen1,m6_aalen_time,by="Var1",all.x=TRUE)
m6_aalen <- m6_aalen[order(m6_aalen$Var2,m6_aalen$Var1),]

m1_aalen$std.err <- sqrt(m1_aalen$total.var)
m1_aalen$lower <- m1_aalen$mean.coef - 1.96*m1_aalen$std.err
m1_aalen$upper <- m1_aalen$mean.coef + 1.96*m1_aalen$std.err
m1_aalen$Var2 <- gsub("`", "", m1_aalen$Var2)
m1_aalen$Var2 <- trimws(m1_aalen$Var2, which = "right")
m1_aalen$Var2.fct <- factor(m1_aalen$Var2,
                            levels=rev(c("Common law",
                                "Population (log)",
                                "GDP capita (log)",
                                "Polyarchy",
                                "Milex per cap",
                                "MID side B",
                                "MID side A",
                                "(Intercept)")),
                            ordered=TRUE)

m2_aalen$std.err <- sqrt(m2_aalen$total.var)
m2_aalen$lower <- m2_aalen$mean.coef - 1.96*m2_aalen$std.err
m2_aalen$upper <- m2_aalen$mean.coef + 1.96*m2_aalen$std.err
m2_aalen$Var2 <- gsub("`", "", m2_aalen$Var2)
m2_aalen$Var2 <- trimws(m2_aalen$Var2, which = "right")
m2_aalen$Var2.fct <- factor(m2_aalen$Var2,
                            levels=rev(c("Common law",
                                "Population (log)",
                                "GDP capita (log)",
                                "Polyarchy",
                                "Milex per cap",
                                "MID side B",
                                "MID side A",
                                "(Intercept)")),
                            ordered=TRUE)


m3_aalen$std.err <- sqrt(m3_aalen$total.var)
m3_aalen$lower <- m3_aalen$mean.coef - 1.96*m3_aalen$std.err
m3_aalen$upper <- m3_aalen$mean.coef + 1.96*m3_aalen$std.err
m3_aalen$Var2 <- gsub("`", "", m3_aalen$Var2)
m3_aalen$Var2 <- trimws(m3_aalen$Var2, which = "right")
m3_aalen$Var2.fct <- factor(m3_aalen$Var2,
                            levels=rev(c("Common law",
                                "Population (log)",
                                "GDP capita (log)",
                                "Polyarchy",
                                "Milex per cap",
                                "MID side B",
                                "MID side A",
                                "(Intercept)")),
                            ordered=TRUE)


m4_aalen$std.err <- sqrt(m4_aalen$total.var)
m4_aalen$lower <- m4_aalen$mean.coef - 1.96*m4_aalen$std.err
m4_aalen$upper <- m4_aalen$mean.coef + 1.96*m4_aalen$std.err
m4_aalen$Var2 <- gsub("`", "", m4_aalen$Var2)
m4_aalen$Var2 <- trimws(m4_aalen$Var2, which = "right")
m4_aalen$Var2.fct <- factor(m4_aalen$Var2,
                            levels=rev(c("Common law",
                                "Population (log)",
                                "GDP capita (log)",
                                "Polyarchy",
                                "Milex per cap",
                                "MID side B",
                                "MID side A",
                                "(Intercept)")),
                            ordered=TRUE)


m5_aalen$std.err <- sqrt(m5_aalen$total.var)
m5_aalen$lower <- m5_aalen$mean.coef - 1.96*m5_aalen$std.err
m5_aalen$upper <- m5_aalen$mean.coef + 1.96*m5_aalen$std.err
m5_aalen$Var2 <- gsub("`", "", m5_aalen$Var2)
m5_aalen$Var2 <- trimws(m5_aalen$Var2, which = "right")
m5_aalen$Var2.fct <- factor(m5_aalen$Var2,
                            levels=rev(c("Common law",
                                "Population (log)",
                                "GDP capita (log)",
                                "Polyarchy",
                                "Milex per cap",
                                "MID side B",
                                "MID side A",
                                "(Intercept)")),
                            ordered=TRUE)


m6_aalen$std.err <- sqrt(m6_aalen$total.var)
m6_aalen$lower <- m6_aalen$mean.coef - 1.96*m6_aalen$std.err
m6_aalen$upper <- m6_aalen$mean.coef + 1.96*m6_aalen$std.err
m6_aalen$Var2 <- gsub("`", "", m6_aalen$Var2)
m6_aalen$Var2 <- trimws(m6_aalen$Var2, which = "right")
m6_aalen$Var2.fct <- factor(m6_aalen$Var2,
                            levels=rev(c("Common law",
                                "Population (log)",
                                "GDP capita (log)",
                                "Polyarchy",
                                "Milex per cap",
                                "MID side B",
                                "MID side A",
                                "(Intercept)")),
                            ordered=TRUE)



m1_plot <- ggplot(data = m1_aalen, aes(x = Time)) +
    geom_stepribbon(aes(ymin = lower, ymax = upper), 
                    fill = "blue", alpha = 0.3) +
    geom_step(aes(y = mean.coef), direction = "hv") +
    facet_wrap(~ Var2.fct, ncol = 4, scales = "free_y",strip.position = "left") +
    theme_bw(base_size = 12) +  # Sets base font size
    theme(strip.text = element_text(size = rel(0.9),
                                    margin = margin(t = 1, r = 1, b = 1, l = 1)),   # Relative to base size
          axis.title.y = element_blank(),
          plot.title = element_text(size = 14, face = "bold", hjust = 0.5)) +
    labs(title = "Overall from 1991")


m2_plot <- ggplot(data = m2_aalen, aes(x = Time)) +
    geom_stepribbon(aes(ymin = lower, ymax = upper), 
                    fill = "blue", alpha = 0.3) +
    geom_step(aes(y = mean.coef), direction = "hv") +
    facet_wrap(~ Var2.fct, ncol = 4, scales = "free_y",strip.position = "left") +
    theme_bw(base_size = 12) +  # Sets base font size
    theme(strip.text = element_text(size = rel(0.9),
                                    margin = margin(t = 1, r = 1, b = 1, l = 1)),   # Relative to base size
          axis.title.y = element_blank(),
          plot.title = element_text(size = 14, face = "bold", hjust = 0.5)) +
    labs(title = "Overall from 1945")


m3_plot <- ggplot(data = m3_aalen, aes(x = Time)) +
    geom_stepribbon(aes(ymin = lower, ymax = upper), 
                    fill = "blue", alpha = 0.3) +
    geom_step(aes(y = mean.coef), direction = "hv") +
    facet_wrap(~ Var2.fct, ncol = 4, scales = "free_y",strip.position = "left") +
    theme_bw(base_size = 12) +  # Sets base font size
    theme(strip.text = element_text(size = rel(0.9),
                                    margin = margin(t = 1, r = 1, b = 1, l = 1)),   # Relative to base size
          axis.title.y = element_blank(),
          plot.title = element_text(size = 14, face = "bold", hjust = 0.5)) +
    labs(title = "Force from 1991")

m4_plot <- ggplot(data = m4_aalen, aes(x = Time)) +
    geom_stepribbon(aes(ymin = lower, ymax = upper), 
                    fill = "blue", alpha = 0.3) +
    geom_step(aes(y = mean.coef), direction = "hv") +
    facet_wrap(~ Var2.fct, ncol = 4, scales = "free_y",strip.position = "left") +
    theme_bw(base_size = 12) +  # Sets base font size
    theme(strip.text = element_text(size = rel(0.9),
                                    margin = margin(t = 1, r = 1, b = 1, l = 1)),   # Relative to base size
          axis.title.y = element_blank(),
          plot.title = element_text(size = 14, face = "bold", hjust = 0.5)) +
    labs(title = "Force from 1945")


m5_plot <- ggplot(data = m5_aalen, aes(x = Time)) +
    geom_stepribbon(aes(ymin = lower, ymax = upper), 
                    fill = "blue", alpha = 0.3) +
    geom_step(aes(y = mean.coef), direction = "hv") +
    facet_wrap(~ Var2.fct, ncol = 4, scales = "free_y",strip.position = "left") +
    theme_bw(base_size = 12) +  # Sets base font size
    theme(strip.text = element_text(size = rel(0.9),
                                    margin = margin(t = 1, r = 1, b = 1, l = 1)),   # Relative to base size
          axis.title.y = element_blank(),
          plot.title = element_text(size = 14, face = "bold", hjust = 0.5)) +
    labs(title = "No force from 1991")

m6_plot <- ggplot(data = m6_aalen, aes(x = Time)) +
    geom_stepribbon(aes(ymin = lower, ymax = upper), 
                    fill = "blue", alpha = 0.3) +
    geom_step(aes(y = mean.coef), direction = "hv") +
    facet_wrap(~ Var2.fct, ncol = 4, scales = "free_y",strip.position = "left") +
    theme_bw(base_size = 12) +  # Sets base font size
    theme(strip.text = element_text(size = rel(0.9),
                                    margin = margin(t = 1, r = 1, b = 1, l = 1)),   # Relative to base size
          axis.title.y = element_blank(),
          plot.title = element_text(size = 14, face = "bold", hjust = 0.5)) +
    labs(title = "No force from 1945")

## knitr AALEN_FIGURE
aalen.plot.1991 <- ggarrange(m1_plot,m3_plot,m5_plot,
#          m2_plot,m4_plot,m6_plot,
#          nrow=2,
          ncol=1 #common.legend = TRUE, legend="bottom2
           )


aalen.plot.1945 <- ggarrange(m2_plot,m4_plot,m6_plot,
                             ncol=1)

ggsave("gc-jg-project/kampala/paper/figures/figure6.png", plot = aalen.plot.1991,
    width = 14, height = 16, device = "png")


ggsave("gc-jg-project/kampala/paper/figures/app-figure-a6.png", plot = aalen.plot.1945,
    width = 14, height = 16, device = "png")

